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Aided Evaluation in Rank-1 LoS Rician Fading 
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Abstract 

We study zero-forcing detection (ZF) for multiple input/multiple output (MIMO) spatial multiplex¬ 
ing under transmit-correlated Rician fading for an Wr x Nj channel matrix with rank-1 line-of-sight 
(LoS) component. By using matrix transformations and multivariate statistics, our exact analysis yields 
the signal-to-noise ratio moment generating function (m.g.f.) as an infinite series of gamma distribution 
m.g.f.’s and analogous series for ZF performance measures, e.g., outage probability and ergodic capacity. 
However, their numerical convergence is inherently problematic with increasing Rician iC-factor, Ar, 
and At. We circumvent this limitation as follows. First, we derive differential equations satisfied by 
the performance measures with a novel automated approach employing a computer-algebra tool which 
implements Grobner basis computation and creative telescoping. These differential equations are then 
solved with the holonomic gradient method (HGM) from initial conditions computed with the infinite 
series. We demonstrate that HGM yields more reliable performance evaluation than by infinite series 
alone and more expeditious than by simulation, for realistic values of K, and even for Ar and At 
relevant to large MIMO systems. We envision extending the proposed approaches for exact analysis and 
reliable evaluation to more general Rician fading and other transceiver methods. 
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I. Introduction 

A. Background, Motivation, and Scope 

The performance of multiple input/multiple output (MIMO) wireless communications systems 
has remained under research focus as the multiantenna architectures that attempt to harvest 
MIMO gains have continued to evolve, e.g., from single-user MIMO, to multi-user and distributed 
MIMO, and, most recently, to massive or large MIMO [[T|, Q, Q, 0, Q, [[^, [|7|. 

As the numbers of transmitting and receiving antennas, herein denoted with Nj and Ar, respec¬ 
tively, have increased in seeking higher array, diversity, and multiplexing gains [[^ pp. 72, 64, 385], 
transceiver processing complexity has also increased. For spatial multiplexing transmission, 
linear detection methods Q, Q, [j^, e.g., zero-forcing detection (ZF) and minimum mean- 
squared-error detection (MMSE), are attractive because of their relatively-low complexity order 
0 (ArAt + ArA| + A|) 0 and their good performance for Ar S> Aj, as the columns of the 
Ar X At channel matrix H tend to become independent [[^. 

For increased practical relevance, MIMO channel model complexity has also been growing, 
and, with it, the difficulties of MIMO performance analysis and numerical evaluation. Thus, early 
ZF research assumed zero-mean, i.e., Rayleigh fading, for the elements of H, which enabled 
relatively simple analysis and evaluation [^, [|^, p0| . Recently, various cases of nonzero-mean 
H, i.e., Rician fading, have rendered increasingly more difficult the analysis and evaluation for 
several transceiver methods [[TT|, [[T2|, [[l^, [[T4|, (HI, |[T6|, |[T7|, (181 , 


i, m- 

Rician fading can occur due to line-of-sight (LoS) propagation, in indoor, urban, and suburban 


scenarios, as shown by the WINNER II channel measurements |22 Section 2.3]. WINNER II 


122 Table 5.5] has also characterized as lognormal the distributions of 1) the Rician A-factor, 
which determines the strength of the channel mean vs. standard deviation (TJ p. 37], and 2) 


the azimuth spread (AS), which determines the antenna correlation [23 p. 136]. An ability to 
evaluate MIMO performance over the range of realistic values of K and AS is useful, e.g., in 


averaging over their distributions, which has only rarely been attempted before [18|. 

Consequently, we focus herein on evaluating MIMO ZF under transmit-correlated Rician 


fading. For tractable analysis we assume as in [ 151, [H| that the EoS or deterministic component 









3 


of H satisfies rank(Hd) = r = 1. Whereas for LoS propagation r ean take any value from 1 

1, small antenna apertures, relatively-low earrier frequeney, or large 


to Nr [24|, [251, 


transmitter-reeeiver distanee, as in eonventional point-to-point deployments [[T|, [ [T^ , are likely 
to yield as outer product of array response vectors [|T| Eq. (7.29), p. 299], i.e., r = 1. 

Our future work shall consider Rician fading with r > 1 for ZF and MMSE. Higher r, which 
improves H conditioning, i.e., MIMO performance, is becoming increasingly more relevant due 
to envisioned EoS millimeter-wave applications [ |24| . MMSE is appealing because it outperforms 
ZE. Also, we shall tackle more general statistical fading models that can characterize more 


modern MIMO deployment types [17|. Einally, for millimeter waves and massive MIMO, we 


shall pursue beamspace channel matrix representation and signal processing [27|, [251, [28 j. 


B. Limitations of Relevant Previous Work on MIMO ZF 

Historically, the study of MIMO ZF commenced with that for uncorrelated Rayleigh fading 
from [[^. The case of transmit-correlated Rayleigh fading was elucidated in Q, [ lOj. For Rician 
fading, previous studies assumed certain values for r and/or proceeded by approximation: 

• Rician fading only for 1) the intended stream, i.e., Rician-Rayleigh fading, which is a 
special case with r = 1, or 2) the interfering streams, i.e., Rayleigh-Rician fading, whereby 
r = At — 1; these cases may arise in heterogeneous networks. Then, we derived in [ fT^ 
exact infinite-series expressions for performance measures, e.g., the average error probability, 
outage probability, and ergodic capacity (i.e., rate [ fT5| , [ [T7| ) — more details below. 

• Rician fading for all streams, i.e.., full-Rician fading, for the special case with r = 1. Early 
works — see [fT8|, [[29| and references therein — used an approximation of the cumbersome 


noncentral-Wishart distribution of with a central-Wishart distribution of equal mean, 
which has yielded a simple gamma distribution for the signal-to-noise ratio (SNR). Then, 
[18 1 , [21 1 reveal that r = 1 does not ensure consistent approximation accuracjQand r > 1 
can render it useless. Recently, bounding techniques have yielded — only for uncorrelated 
fading — the simple sum rate bounds in [ [T5| Eqs. (55)-(58)] that become accurate at high 
SNR. 


*Only very careful usage in |l8) helped average the performance over WINNER II distrihutions of K and AS for r = 1. 
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Rician fading, Vr = 1, ■ • • , Nj. For this most general ease, exaet sum-rate expressions for 
A^r —)■ oo and approximations for finite A^r were derived in 


For Rieian-Rayleigh fading, we have reeently analyzed and evaluated ZF exaetly in [19|, 


|20|. In [19|, we expressed the SNR moment generating funetion (m.g.f.) in terms of the 
eonfluent hypergeometrie funetion iFi(-,-,(t) mi Eq. (31)], where a (x KN^^Nj. Thereafter, 


its well-known expansion around cxo = 0 [19, Eq. (30)] yielded an infinite series of gamma 


distribution m.g.f.’s [ 19, Eq. (37)]. Einally, inverse-Eaplaee transformation and integration yielded 
analogous series for the SNR probability density funetion (p.d.f.), average error probability, 
outage probability, and ergodie eapaeity 03 Eqs. (39), (58), (69), (71)]. However, beside 
eomplieating the analysis, the Wishart distribution noneentrality indueed by Rieian fading also 
leads to numerieal divergenee for these series with inereasing K, Nr, and A^x [ [T9l Seetion V.E]. In 
|20|, we overeame this limitation by using the faet that iFi(-, •, a) is a holonomic functiot^ i.e., it 


satisfies a differential equation [20, Eq. (27)] with polynomial eoeffieients with respeet to (w.r.t.) 
a. Starting from this differential equation, a diffieult by-hand derivation produeed differential 
equations for the SNR m.g.f. and then for the SNR p.d.f., via inverse-Eaplaee transform. 
Thereafter, we eomputed reliably the p.d.f. at realistie values of K — but only for relatively small 
A^r and Nj — by numerieally solving its differential equations from initial eonditions eomputed 
with the infinite series for small K. This approaeh is known as the holonomic gradient method 
(HGM) beeause, at eaeh step, the funetion value is updated with the differential gradient 
See. rV.B]. Einally, in [ |20| , the SNR p.d.f. eomputed with HGM was numerieally integrated to 
evaluate performanee measures, i.e., the outage probability and ergodie eapaeity. 

Thus, on the one hand, our exaet studies for r = 1 in [ [T^ , [20[ are limited by the following: 

• Nonfull-Rieian (i.e., only Rieian-Rayleigh) fading assumption. 

• Tedious by-hand derivations of the SNR m.g.f. and p.d.f. differential equations. 

• Time-eonsuming numerieal integration of the p.d.f. for performanee measure evaluation. 

• HGM not tried for large A^r and A^x, e.g., as relevant for large MIMO systems Q. 

On the other hand, only approximations exist for full-Rieian fading and r = 1 [|T5|, [[T7|, [|T^. 


^Other examples: rational functions, logarithm, exponential, sine, special functions (orthogonal polynomials, Bessel 


p. 


41]). 
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C. Problem Tackled in the Current Work; Exact Analysis and Evaluation Approaches 

To the best of our knowledge, the performanee of MIMO ZF has not yet been studied exactly 
under full-Rician fading even for r = 1. We pursue this study herein, as follows. 

First, upon applying a sequenee of matrix transformations and results from multivariate statis- 
ties, we obtain several theoretieal results that help express exactly the SNR m.g.f. as an infinite 
series with terms in iFi(-, •, •). Thus, the m.g.f. can be rewritten as a double-infinite series of 
gamma distribution m.g.f.’s, which readily yields analogous series for the SNR p.d.f. and for 
the performance measures. Then, they are recast as a generic single-infinite series. However, its 
truncation is found to incur numerical divergence with increasing K, TYr, and N^. Consequently, 


as in 1191, it is necessary to derive satisfied differential equations and apply HGM. 


Because by-hand derivation of differential equations for our generic series appears intractable, 
we resort to a novel automated derivation approach using the HolonomicFunct ions package 


written earlier by one of the authors |30|, [311 and implementing recent advances in computer 
algebra. It exploits, for holonomic functions, closure properties [|^ Section IV.C], [ |^ , the 
algebraic concept of Grobner |33|, and creative telescoping algorithms [^, Ch. 3] to 

systematically deduce differential equations for their addition, multiplication, composition, and 
integration. This computer-algebra-aided approach readily yields differential equations not only 
for the SNR m.g.f. and p.d.f., but also for the outage probability and ergodic capacity. 

Finally, we evaluate ZF performance measures by HGM, i.e., by solving the obtained differ¬ 
ential equations starting from initial conditions computed with the infinite series. 


D. Contributions 

Compared to previous MIMO ZF work by us and others, herein we: 

• Tackle full-Rician fading with r = 1 in a new exact analysis that reveals that the SNR 
distribution is an infinite mixture of gamma distributions. This SNR distribution yields 
insight into the effect of channel matrix statistics (mean, correlation) on performance, and 
helps reassess the approximation with the gamma distribution we studied in [18|, pT|. 


^Buchberger’s algorithm 


32 


for Grobner basis computation specializes, for example, to the Euclidean algorithm when applied 


to univariate polynomials, and to Gaussian elimination when applied to linear polynomials in several variables |33| . Grdhner 
bases have helped solve communications optimization problems cast as systems of polynomial equations, e.g., for interference 
alignment |34|, coding gain maximization in space-time coding (35); other relevant applications are listed in ||33). 
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• Use computer algebra to automate deductions of differential equations also for performance 
measures and, thus, also avoid time-consuming numerical integration of the SNR p.d.f.. 

• Demonstrate that HGM yields accurate performance evaluation for realistic values for K, 
and even for large TYr and Nj, unlike the infinite series alone and faster than by simulation. 

• Exactly average the ZF performance over WINNER II distributions of K and AS. 


E. Paper Organization 

Section |l^ introduces our model. Section |l^ employs matrix transformations and multivariate 


statistics to express exactly the m.g.f. of the ZF SNR. Section IV derives a generic infinite 
series for the SNR m.g.f. and p.d.f., as well as for ZF performance measures. Section |V] 
describes the automated derivation of differential equations, which has been accomplished with 
HolonomicFunctions commands as shown in [^. Finally, Section [Vl| presents numerical 
results obtained by simulation, series truncation, and HGM. The Appendix shows some proofs 
and derivation details. 


F. Notation 


• Scalars, vectors, and matrices are represented with lowercase italics, lowercase boldface, and 
uppercase boldface, respectively, e.g., y, h, and H; the statement H = Ar x At indicates 
Ar rows and Aj columns for H; zero vectors and matrices of appropriate dimensions are 
denoted with 0; superscripts and stand for transpose and Hermitian (i.e., complex- 
conjugate) transpose; I at is the A x A identity matrix. 

• [-ji is the Ah element of a vector; [■]* ,, and [-j.j indicate the i,jth element, Ah row, 

and jth column of a matrix; ||H|p = is the squared Frobenius norm. 

• i = 1 : N stands for the enumeration i = l,2,...,A;(8) stands for the Kronecker product 
p7] p. 72] ; oc stands for ‘proportional to’; ^ stands for logical implication. 

• h ~ CAtvr (hd, R) denotes an Ar x 1 complex-valued circularly-symmetric Gaussian vector 

with mean hd and covariance matrix R; an Ar x At complex-valued circularly-symmetric 
Gaussian random matrix with mean Hd, row covariance Iat^, and column covariance Rt, i.e., 
a matrix whose vectorized form is distributed as vec(H'^) CNATjjAfT )5 Itvr ^ Rt)> 

is denoted herein as H ~ CAfArR,ArT (Hd, Iatr 0 Rt), based on the definition from pA| ; 
subscripts -d and t identify, respectively, deterministic and random components; subscript 





7 


■n indicates a normalized variable; E{-} denotes statistieal average; r(iV, Fi) represents the 
gamma distribution with shape parameter N and seale parameter Fi; Xmi^) denotes the 
noneentral chi-square distribution with m degrees of freedom and noneentrality parameter 
5\ Xm denotes the central ehi-square distribution with m degrees of freedom; B(A^, M) rep¬ 
resents the eentral beta distribution with shape parameters N and M; B(iV, M, x) represents 
the noneentral beta distribution with shape parameters N and M, and noneentrality x. 
iFi(-; ■; ■) is the confluent hypergeometric function [38 Eq. (13.2.2)]; (iV)„ is the Poehham- 
mer symbol, i.e., (A^)o = 1 and {N)n = N{N -y 1)... {N + n — >1. 

d^g{t,z) denotes the kth partial derivative w.r.t. t of funetion g{t,z). 


II. Model and Assumptions 
A. Received Signal and Fading Models 

We eonsider an uneoded point-to-point uplink MIMO spatial multiplexing system over a 
frequeney-flat fading ehannel [[^ Chs. 3, 7]. There are Nj >2 and Ar > Ap antenna elements 
at the transmitteij^ and reeeiver, respectively. For the transmit-symbol veetor denoted with 

y = {yi 2/2 ■ ■ ■ UNj)^ = At X 1, (1) 

the stream of eomplex-valued symbols Pi from antenna i is referred to as Stream i. Without loss 
of generality, we eonsider Stream 1 as the intended stream (i.e., whose symbol is deteeted, and 
whose deteetion performance is analyzed and evaluated), and the remaining 

Aj = Aj — 1 (2) 

streams, i.e.. Streams i = 2 : Ap, as interfering streams. The number of degrees of freedom is 

A = Ar — Aj = Ar — Ap -|- 1. (3) 


Then, the reeeived signal veetor ean be represented as 

r = W ^ Hy + n = Ar X 1, 
V -/V7 


(4) 


'’For Ni = 1 and maximal-ratio combining (MRC), we obtained a simple SNR m.g.f. expression for Rician fading in jl9[ 
Eq. (36)]. 



where ^ is the energy transmitted per symbol (i.e., per antenna), and n ~ C!NAr^( 0 , Iat^) is 


the additive noise. Then, the per-symbol transmit SNR is 

y _ E, 1 


(5) 


Finally, we assume that the eomplex-valued ehannel matrix H = TVr x Nt is Gaussian (more 
details follow below), has rank N^, and is perfeetly known at the reeeiveij^ With its deterministie 
and random eomponents denoted as Hd and Hr, respeetively, we ean write 


H = Hd + Hr = 


K 


H 


K + 1 ’ V A' + 1 

where Hd,n and Hr.n are the eomponents of H normalized as 

d,n|p = E{||H,„||2} = ATrATt, i.e., E{||Hf} = N^Nj, 


Hrn, 


IH 


and K, known as the Rieian iF-faetor, is deseribed by 


K = 


iHdiP 


K 

K+1 


IH 


d,n I 


( 6 ) 


(7) 


( 8 ) 


E{||Hr||2} ^E{||Hr,||2}- 

Then, Ff = 0 yields full-Rayleigh fading, i.e., | [H].^. | is Rayleigh distributed Vz, j, as assumed 
in [[^, 0, [10|. Further, the ease when iT 7 ^ 0 and in Hd,n only eolumn [Hd,n].,i is nonzero is 

1. Finally, herein, the ease when iT 7 ^ 0 


referred to as Rieian-Rayleigh fading, as in [19|, 
and eaeh eolumn of Hd,n has at least one nonzero element is referred to as full-Rieian fading. 

We assume that Hd arises due to LoS propagation between transmitter and reeeiver. Then, 
if the transmitter-reeeiver distanee is mueh larger than the antenna interelement spaeing, Hd 
ean be represented as the outer produet of the array response veetors for the reeeiving antenna, 
a = A^r X 1, and transmitting antenna, b = A^t x 1, i-O-, 0 Eq. (7.29), p. 299] 


Hd = ab^ = a ( 6 * 6 * 




(9) 


whieh reveals that Hd has rank r = 1 and eolumns given by hd,* = a 6 *, z = 1 : A^t- 
Remark 1. We may assume that ||a|| = 1 if we scale b according to 

Nj Nj Nj 


|bf = ^ \hi\^ = ^ ||af ^ ||hrf,i| 


2 = 1 


2 = 1 


= 1 


2=1 


2 = ||H,f 

A + 1 


(10) 


^ZF for imperfectly-known H can be studied, e.g., with the effective-SNR approach we described in |l8|. 
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For a tractable analysis, we assume zero row eorrelation (i.e., reeeive-antenna eorrelation) for 
H. On the other hand, we assume, as in 0, p^ , p^ , p0[ , that any row of Hr,n has the same 
distribution C?sfAr^(0, Rt), so that any row of Hr has the same distribution C?sfAr^(0, Rt,j^) with 

1 


R„. = H„} = ^ ^ 


Rt- 


( 11 ) 


1 /2 

Thus, we ean write Hr = HwRjy^ with H„ 

C^AfR,AfT (Hd, IaTr ® Rt,x)- 

Matrix Rt is determined by antenna interelement spaeing and AS, i.e., the ‘standard deviation’ 






( 0 , Inr ( 8 ) lN^), so that H = Hd + Hr 


of the power azimuth spectrum [23 p. 136]. When the latter is modeled as Laplaeian, as 
reeommended by WINNER II Rt ean be eomputed from the AS with [|23| Eqs. (4-3)-(4-5)]. 


Remark 2. WINNER II modeled the measured AS (in degrees) and K (in dB) as random variables 


with scenario-dependent lognormal distributions [22 Table 5.5] [18 Table 1]. Thus, herein, we 


attempt to evaluate ZF performance for AS and K values relevant to these distributions. 


B. Matrix Partitioning Used in Analysis 

To study Stream-1 deteetion performanee, we shall employ the partitioning 

H=(hi H2) = (hd,l Hd,2) + (hr,l Hr, 2 ), 


( 12 ) 


where hi, hd,i, and hr,i are Ar x 1 veetors, whereas H 2 , Hd, 2 , and Hr, 2 are Nr x Ni matriees. 
We shall also employ the corresponding partitioning of the eolumn eovariance matrix: 


= I I = I ^’-21 I (13) 

Rr.irai R'T,X22 

Remark 3. Herein, we consider full-Rician fading with r = rankifUd) = ranfc(H^, 2 ) = 1, 
whereas in we considered its special case of Rician-Rayleigh fading, i.e., rankifAf) = 

1, but rank(Ylii^ 2 ) = 0. Thus, the results obtained herein specialize to those in 0 0 when 
we reduce to 0 the vector formed with the last N/ = Nt — 1 elements ofh, i.e., the vector 


b — (62 ... 


(14) 
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III. Exact Analysis of ZF SNR 
A. ZF SNR as Hermitian Form 


Given H, ZF for the signal from Q refers to symbol detection based on the operation 


n 


^ yr; L J yivo 


(15) 


Based on ( [15] ) and [ [T^ , p^ , the SNR for Stream 1 can be written as the Hermitian form below: 

r. 


7i = 


[(HHH) 


= r.h" hi, 


(16) 


Ji,i 


=Q2 


where Q 2 = Nr x Nr is idempotent and of rank N. 


Remark 4. The following transformations do not change the ZF SNR in {16): 

• Row transformations o/H with unitary matrices, because they do not change 

• Column transformations 0 /H 2 with nonsingular matrices, because they do not change Q 2 . 
Several such transformations, shown below, help derive the exact SNR distribution. 


B. Row Transformation F = VH That Zeroes Rows [F^]j„ i = 2 : Nr 


If we make the substitution H = V'^F, with unitary V = Nr x Nr, in (16) and partition 


(fr, 


1 r 1-^27 


(17) 


(18) 


according to (12) the matrix 

F = \U = NrxNt 

= (fi F2) = (fd,l Fd,2) 

the ZF SNR Hermitian form in ( [T^ becomes 

7i = Fsff Q2fi, 

with 

Q2 = 1n,-F2{F^F,)-^FI (19) 

Choosing the first row of the unitary matrix V as [V]i, = we conveniently obtain 
[Fd]l,.^([V]l,.a)bH = ||a||2b^^ = b^ 

[Fd]*,.i ([V],,.a)b^'=0, z = 2:Nr, 


=0 
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i.e., fdj (6* 0 


0 )"'', j — 1 : -^T- 


( 20 ) 


Theorem 1. The m.g.f. of the SNR conditioned on Q 2 can be written, simply, as 

1 


^7iIQ2(s) = E^i{e"^'IQ2} = 


{l-TisY 

with scalar Fi and function /i(s) defined in the proof below. 


exp{/i(s)[Q 2 ]i,i}, 


( 21 ) 


Proof: 

Because the column covariance of F = VH is the same as that of H, i.e., partitioned 


as in ( |13[ ), and because fi = A^r x 1 and F 2 = A^r x Ni from the partitioning of F in ( |17[ ) are 
jointly Gaussian, the distribution of fi given F 2 is given by Appendix], [19, Eqs. (12)-(16)] 


filFs 


CK 


atrI (fd,i - Fd,2r2,i)^ +F2r2,i, ([Rt,k]ij) > 

=|J.=Nf^xl 


( 22 ) 


with 


r2,i - - ivi 

1 1) = ''"T.ATii — ’'T.Kai ^,^"22 


-1 


= At X 1, 


Ko 


(23) 

(24) 


Then, it can be shown by substituting (|22j) into (181 and further manipulating as in [ lOj, [ 19|, 


that the SNR conditioned on Q 2 from ( |T^ can be written as the Hermitian form 

71 IQ 2 


Fi 

fi 


= rif|^Q 2 fi, with 

r, 

Ay,,’ 

~ CAtvr 




Nk ) , 


122 [ 


J20l 


(25) 

(26) 

(27) 

(28) 


II - fd,i - Fd,2r2,i - (6] - b^r2,i 0 ... 0)"^ = (/n 0 ... 0)"^, 

i.e., row transformation F = VH yielded a single nonzero-mean element in fi, which simplifies 
the ensuing analysis. 


The Hermitian form in fi from (25) helps cast the m.g.f. of the SNR given Q 2 as [ 19 Eq. (20)] 

exp { - xiu^ [Ijv^ - - risQ 2 )“^] u} 


M 


71 IQ: 


(^) = 


det (Iatj, - FisQs 


with 




(29) 


(30) 
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^=^ = (1 0 


Imr — (Iatr — FisQ 


0 )T 

l-Fis 


Q 2 


(31) 

(32) 


Above, (32) follows by using the eigendeeomposition of Q 2 . The desired m.g.f. expression 


in (21) follows by substituting (32) into (29) and defining /i(s) = 


C. Partial Column Transformations That Help Rewrite [Q 2 ]i,i Conveniently 

1) Unitary Transformation E 2 = F 2 V That Zeroes Elements J = 2 : Np Making 

the substitution F 2 = F 2 V'^, with unitary V = iVi x Ni, in (19) yields 

Q2 = IiVR-E2(EHF2)-'F2^ (33) 

Based on ( [TT] ), we ean write 

E 2 = F 2 V = Fd^2^ + Fr^2^ = E(1_2 + Er^2 = -^R ^ -^I- (34) 

Setting [V],_i = b/||b|| simplifies the ensuing SNR analysis as it zeroes [Ed, 2 ]ij, j = 2 : Np 


Ed,2 = Ed, 2 V ^ 



^ [V].,2 ---[VU 


= lib 



(35) 


2) Nonsingular Transformation That Decorrelates the Columns of E 2 : For the eolumn eor- 


relation of E 1 . 2 from (34), i.e., for 


Te{E,“2E,,,} = Te{(F,2V)"(F,,2V)}®V"Rt,k„V, 

JVr JVr 


let US eonsider the Cholesky deeomposition ||^ See. 5.6] 

VHRt,x,,V = AA^ 


(36) 


(37) 


where A = Aj x Aj is upper triangular with real-valued and positive diagonal elements. 

Then, eonsidering matrix Evi ,,2 ~ C3\fArR,Ar, (0 , IaTr ® Iatj), we ean write ([^ based on ([^ 


and (36) as 


E 2 — Ed ,2 + Ew, 2A'^ — (Ed,2A + Ew,2) A'^. (38) 

Thus, by transforming the eolumns of E 2 with A^'^, we obtain 

G 2 = E 2 A = Ed,2A + Ew^2 = Ar X Ai, 


(39) 
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whose mean can be written, based on (35) and the fact that A is lower triangular, as 

1 0 


Gd,2 = Ed,2A-^^ = ||b||[A-^^]i,i 


0 0 


(40) 


Using (37), the properties of A, and the choice [V],^i = b/||b||, the squared norm of Gd ,2 can 
be written as 

:C2=l|G42||"=b"ltF3f^^b. (41) 

Remark 5. For Rician-Rayleigh fading, Remark^revealed that b = 0, which by {41) implies 
X 2 = 0. On the other hand, for full-Rayleigh fading, implies that also Xi = 0. 


Thus, column transformation ( [39] ) yielded G 2 with uncorrelated columns and mean given by 

, (i.e., real-valued) for i = j = 1, 


[G, 




(42) 


0 , otherwise. 


Substituting E 2 = G 2 A'^ into (33) yields 


Q2=lAr,-G2(GHG2)-'G 


ir^H 
2 • 


(43) 


The simple statistics of G 2 (vs. F 2 ) help simplify our SNR distribution analysis, as shown 
below. 


3) QR Decomposition: Finally, by substituting in ( [43] ) the QR decomposition [j^ Sec. 5.7] 

G 2 = U 2 T 2 , (44) 


where U 2 = Ar x Ni satisfies U 2 U 2 = Iatj, and T 2 = Aj x Aj is upper triangular with 
real-valued and positive diagonal elements, we can write Q 2 simply as 

Q 2 = In, - U 2 T 2 (THT 2 )-'THuH = - U 2 U 2 ^ (45) 

This helps write [Q 2 ]i,i for the m.g.f. in ( jlTj ) solely in terms of the first row of U 2 as 

[Q2 ]i,1 = 1 - [U2]l,. ([U2]l,.)'' = 1 - (l[U2]l,ir + l[U2]l,2r + ' ' ' + l[U2]l,iV.r) 

= (i-i[U2]i,ir) 

"-V-' 

=/3i 

=/32 


|[U2]l,2p + --- + |[U2]l, 


Ni\ 


l[U 


2 1,1 


(46) 
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D. Principal Analysis Result: Exact M.G.F. Expression of the Unconditioned SNR 


The above transformations have helped write the eonditioned-SNR m.g.f. from (21) as 

1 


ATyi (s I /?!, 132 ) — 


(1 - Fis) 


N 


exp{fi{s)/3i/32}. 


(47) 


In order to express the unconditioned-SNR m.g.f., we need to average (47) over the distributions 
of (3i and (32, whieh are elueidated in the following two lemmas. 


Lemma 1. Random variable (3i from ( |46P is distributed as 

( 3 i ~ B{Nr - 1 , 1 , 0 : 2 ). 

Proof: See Appendix 


(48) 


Lemma 2. Random variable (32 from (46) is distributed as 

/32~B(iV,iVz-l), 


i.e., has m.g.f [19 Eq. (30)] 


= iFi{N;Nr-1;s), 


(49) 


(50) 


and is independent of (3i. 

Proof: See Appendix 

Theorem 2. The m.g.f. of the unconditioned ZF SNR under full-Rician fading with r = 1 is 


Xi,X2) = 


E 


6 ^2 


> 12 ! 


iFi \N]n2 + Nr] 




-Xi 


(51) 


1 — Tis 

Proof: Due to limited space, we only outline the proof: it follows by successively averaging 


the m.g.f. of the conditioned SNR in (47) over the distributions of the independent random 


variables (3i and (32, and by exploiting (50), (53), and (83). 


E. Effects of Channel Matrix Statistics on SNR Statistics 

For Rician-Rayleigh fading (i.e., for X 2 = 0), the SNR m.g.f. from ( [ST] ) reduces to [19 
Eq. (31)] 

1 


M.^^(s;xi) = 


(l-Fis) 


N 




(52) 
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Then, for 71 , the first two moments, varianee ¥{ 71 } = IE{ 7 i} — (E{ 7 i})^, and amount of fading 


■^{ 71 } = ¥{ 71 }/ (E{ 7 i}) , i.e., SNR statisties, have been expressed in [19 Table I]. Beeause 


the SNR m.g.f. for full-Rician fading from (51) is a weighted infinite series of SNR m.g.f.’s 


for Rieian-Rayleigh fading from (52) with TYr replaeed with TYr + ^ 2 , expressing E{ 7 i} and 
E{ 7 ^} for the former from those for the latter from [19, Table I] is trivial. Expressing ¥{ 71 } 


and A{ 7 i} based on ( [sT] ) and [19 Table I] is not trivial. 

On the one hand, the effect of X 2 on SNR statistics is not readily discernible from (51 ) and [ [T^ 
Table I]. On the other hand, ( [sT] ) and [ 19, Table I] reveal that Ejyi} increases with N from ([^, 
Ti from (26), and xi from ([^. Further, note that it can be shown that xi oc ||hd,i — Hd, 2 r 2 ,i||- 


Thus, the performance of ZF for full-Rician fading with r = 1 is worst when the channel matrix 
statistics satisfy condition hd^ = Hd, 2 r 2 ,i, and it improves with increasing ||hd 7 — Hd, 2 r 2 ,i||- In 
0 . where we studied full-Rician fading irrespective of r, we had noticed (e.g., by comparing 


[21 Figs. 1, 2]) that ZF performed worst for hdj = Hd, 2 r 27 . 


Remark 6. Note that Va; 2 , if hj,i = 2 ^ 2 , 1 , i-S-, Xi = 0, then the m.g.f. in ( [57[ ) reduces to 

the gamma m.g.f. = (1 — Tis)"^. On the other hand, the gamma distribution with 


m.g.f. M{s) = (1 — sTi) ^ and Ti obtained as in (26) from IIt,k = ^t,k + has 


previously been employed to approximate the actual ZF SNR distribution for Rician fading, 


irrespective ofr — see [18], [29] and references therein. Interestingly, condition h^i = 2 ^ 2,1 


yields Ti = Ti, rendering the approximation exact — see [21 Corollary 4]. 


The above have yielded the following insights. 

Remark 7. For ZF under full-Rician fading with r = 1, condition h^i = Hj 2 i' 2 ,i yields: 1) worst 
performance; 2) full accuracy for the gamma distribution previously employed to approximate 
the SNR distribution. 


IV. Exact Infinite Series Expressions eor ZF Performance Measures 
A. Infinite Series Expansion of iFi(-; ■;(j) Around cxo = 0 

Using the well-known infinite series expansion around cxo = 0 (13 Eq. (30)] 


iFi(iV;iVR;cr) = 5 ^ 




(53) 
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the SNR m.g.f. from ( p^ for Rician-Rayleigh fading can also be written as [19, Eq. (37)] 




E 


rii 


- 1 ) 


mi 


(iVR)n, ni! W 

m =0 ^ ^ m^=0^ ^ V > 

’V" 

= ]\4n^ ,m^ (- 5 ) 

where mi(s) is the m.g.f. of a random variable distributed as r(iV + rii — m, Ei^ 


(54) 


Theoretically, (53) converges Va. Nevertheless, the computation of (53) by truncation incurs 
inherent numerical convergence difficulties with increasing a 03- Consequently, the computation 
of ensuing measures, e.g., the ZF SNR p.d.f., becomes nontrivial at realistic values of K, as 


revealed in [19|, [20|. Similar difficulties arise also for the case studied herein, i.e., full-Rician 


fading with r = 1 , upon infinite series expansion of iFi(-;-;(t) in the SNR m.g.f. expression 
from dSTj ), as discussed below. 


B. Exact Double-Infinite Series for M.G.F., P.D.F, and Performance Measures 

By substituting (^) into ([ST]) and proceeding as for ([54|), the SNR m.g.f. becomes 


Xi,X2) = ^ 


{NU xT xT 


ni 

(^^2 + iVR)ni ni ! 7 ^ 2 ! 

ni=0 772=0 mi=0 ^ 


E 




(55) 


Using the m.g.f.-p.d.f. Laplace-transform pair corresponding to r(iV + rii — m, Fi), i.e., 

,mi (t’) 


Pni ,mi (f ) 


(1 - sri)^+”i-”^i ’ 


[(Ar + ni-mi)-l]!rf+"i-”^i’ 

the ZF SNR p.d.f. corresponding to (|5^ can be written, analogously, a^ 




V xT I "1 




—^ (no + Nk)„, nf. no\ ^ \mi 

ni= 0 n 2=0 ^ 1 ^ mi=0 ^ ^ 


(0 • 


(56) 

(57) 

(58) 


■ V' 

=p«l (t) 


By integrating (58), the Stream-1 outage probability at threshold SNR r and the ergodic 
capacity (i.e., rate) are exactly characterized by analogous infinite series, i.e.. 


Po(xi,X2) = / xi, X 2 ) dt 

Jo 


®An alternate p.d.f. expression, for real-valued H, appears 


(59) 


39 


Eqs. (18), (31)]. 
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= 


EE 


(N) 


ni 


m n2 

Xi X 2 


ni 

~^n ('^2 + ^R)ni ni! Tla! V"^l 

m=0n2=0 mi=0 ^ 


E 


iWi p 

■^J o,ni ,mi ? 


= ^^o,ni 


C(a:i,X 2 ) = / ln(l+ t)pT,i(t;a;i,a; 2 )dt 


= e-"^ 


^0 

00 00 

EE 

("2 + "i- " 2 ! ^0 Wi 


(AT) 


ni 


ni no ^1 

Xi 0:2 


E 


/_1 \mi^ 

V -*-/ ^ni ,mi 5 


whereQ 


=c„i 


|(Ar + ,ii- ,ni)-l]! ’ 

1 /■“ 

1 ^ / ln(l -|- ^)Pni,mi(^) dt. 

In 2 Jo 


(60) 

(61) 

(62) 


(63) 

(64) 


Finally, the approach in Section V.A] can help express also the average error probability as 
an infinite series analogous to ( [^ and ^2) . 

On the other hand, the previously employed approximating gamma distribution for the ZF 
SNR mentioned in Remarkyields simple performance measures expressions similar to (63) 


and (64). 


C. Generic Single-Infinite Series for M.G.F., P.D.F., and Performance Measures 


Because ( [55| ), ( [58[ ), ( | 6 ()| ), and ( |62| ) are analogous, we may represent them as the generic double 
infinite series 


00 00 


Hxux,) = e-“ y; y; 


“(i n (^R + ^2)ni nf. n 2 \ 

ni=0 n2=0 ^ ^ ^ 


ni 5 


(65) 


where Hn^ stands for Mnfis) from (55), Pn^if) from (58), Po,ni from (60), and Cn^ from (62) 


Thus, the dependence of h{xi,X 2 ) on s for the m.g.f. or t for the p.d.f. is not explicitly shown 


in ( [65] ), for simplicity. 

Numerical results not shown due to length limitations have revealed that increasing K, Nr, 


and Nj yield increasingly problematic numerical convergence for series (65). This is explained 


by: 1 ) the fact that (|55|) has been obtained from (|51|) by replacing 1 Fi ( A^; n 2 + A^r; 


l-Fis" 


^ 'y{k,x) = Jq F * dt is the incomplete gamma function jssj p. 174]. Integral i64i is expressed in Eq. (73)] 
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with its expansion around xi = 0 from ( [5^ ; 2) the faet that xi is inereasing beeause of the 
following proportionality, proved in Appendix 

Xi oc KN^Nj. (66) 

Appendix 1^ also shows that X 2 oc KN^Nt. Then, the expressions for xi and X 2 dedueed there 


in (95) and (96) ean be used to show that their ratio ci = ^ is real-valued, positive, and 

_ _ 


independent of K and Ar. Finally, unshown numerieal results have revealed that Ci oc l/Nj. 


These eonsiderations suggest substituting X 2 = z and xi = ciz in the generie series in ( [65] ), 
whieh yields the following result. 


Lemma 3. For X 2 = z and Xi = Ciz, series ( |65| ) can be recast as the single-infinite series 

oo n 


n\ Wm „ 

' / ^ \ J / AT I ^ ^rnCi . • 

—^\m/ Ao + n — m)m n\ 

n=0m=0 \ ^ ^ 


(67) 


= Gn 


Derivatives of h{z), required below for HGM, are given by 


^ fk 


1=0 




oo 

n=l 


^n—l 


{n — l)\ 


( 68 ) 


Proof: The proof of the first part is not shown, due to simplieity and length limitations. 
The seeond part follows from ([67]) based on Leibniz’s formula [[38[ Eq. (1.4.12), p. 5]. ■ 


Numerieal results shown later reveal that the truneation of (67) still does not eonverge numeri- 
eally for praetieally relevant values of K, Ar, and Aj. Therefore, we shall endeavor to eompute 


it by HGM, as done for Rieian-Rayleigh fading in [20j to eompute the SNR p.d.f. series 


dedueed from (54). Reeall that HGM evaluates a funetion at given values for its variables by 
numerieally solving its differential equations starting from initial conditions, i.e., known values 


of the function and required derivatives, at another point [20 Sec. IV.B]. Thus, HGM requires 


differential equations. Note that, making the substitutions X 2 = z and xi = ciz and regarding 


Cl as a constant factor, conveniently reduces the number of variables in generic series (67). For 
example, when cast for the m.g.f., the series is only a function of s and z. 

Differential equations were derived by hand, with difficulty, for the ZF SNR m.g.f. and p.d.f. in 


[20 Eqs. (32), (42)] for the Rieian-Rayleigh fading case, based on the SNR m.g.f. expression 


shown here in (52) and the differential equation satisfied by iFi(A; Ar; a), i.e., [20, Eq. (27)] 
a - iFf^(A;AR;a) + (AR-a)- iFW(A; Ar; a) - A ■ iFi(A; Ar; a) = 0. 


(69) 


















19 


For the full-Rician fading case with r = 1 studied herein, the new SNR m.g.f. expression 
in ( [ST] ) comprises an extra sum compared to (52). On the other hand, (67) yields the following 
complicated SNR m.g.f. expression: 


z) = e 


M.n{s)= 

mi=0 


EE 

n=0 m=0 

m 
rrii 


/ 

t (N)raMm(s)c^ Z^ 

(70) 

Vmy 

' (Ar + n- m)m n\ ’ 

(-1) 

mi ^ 

(1 - ■ 

(71) 


Because the by-hand derivation of differential equations w.r.t. s and 2 ; satisfied by 
described by a or ( |70l ) is not tractable, we shall apply instead the automated approach 


described below, based on the generic expression (67). The derivation of differential equations 


satisfied by Po{z), and C{z) can be automated as well, based on: 1) their generic 


expression (67); or 2) the Laplace-transform relationship between M^^{s-,z) and p^^{t-,z), and 
the integral relationships of z) with Po{z) and C{z). We shall employ the latter approach 
because it is more general. 


V. Computer-Algebra-Aided Derivation oe Differential Equations for HGM 
A. Holonomic Functions, Annihilator, Grobner Basis, and Creative Telescoping 

A function is holonomic w.r.t. a set of continuous variables if it satisfies for each of them a 
linear differential equation with polynomial coefficients. A function is holonomic w.r.t. to a set of 


discrete variables if the associated generating function is holonomic in the previous sense [|^ 
Sec. IV.C] 


p. 17]. For example, iFi(A^; Ar; a) is holonomic w.r.t. a because it satisfies 
differential equatiorp] (|6^. In other words, iFi(A;Ar;ct) is annihilated by the differential 
operator ad1 + (Ar — a)da — A. The (infinite) set of all operators that annihilate a given 
holonomic function is called its annihilator [ |^ p. 18]. 

Holonomic functions are closed under addition, multiplication, certain substitutions, and taking 
sums and integrals [20|, |3^. Consequently, functions z), z), Po(z), and C(z), cast 


as in (67), are holonomic. The fact that the closure properties for holonomic functions can 
be executed algorithmically provides a systematic way of deriving the differential equations 
required for HGM, by starting with the annihilating operators of the comprised “elementary” 


*Note that iFi(N; Nr; a) is also holonomic w.r.t. N and Nr. 
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holonomic functions in (67). A key ingredient for algorithmieally exeeuting elosure properties is 


the algebraie eoneept of Grobner basis, whieh provides a eanonieal and finite representation of 
an annihilator and helps deeide whether an operator is in an annihilator. For details on Grobner 


bases theory, eomputation, and applieations see [32|, [30|, 


, [33| and references therein. 


While many holonomie elosure properties require, basieally, only linear algebra, eomputing the 
annihilator for a sum or integral of a holonomie funetion is a more involved task. For example, 
one ean employ the creative telescoping teehnique: given an integral F{x) = f{x,y) dy, 
ereative teleseoping algorithmieally finds in the annihilator of f{x,y) a differential operator of 


the form P{x,dx) + dy ■ Q{x,y,dx,dy). Then, using the fundamental theorem of calculus [38 
p. 6] and differentiating under the integral sign reveal^ P{x, d^) as an annihilating operator 


for F{x) [30 p. 46]. Several ereative teleseoping algorithms are deseribed in [30 Ch. 3] 


B. The HolonomicFunctions Computer-Algebra Package 

This freely-available eomputer-algebra paekage, written earlier in Mathematica by one 


of the authors, is deseribed, with numerous examples, in [31 j. Its eommands implement: 1) the 
eomputation of Grobner bases in operator algebras, 2) elosure properties for holonomie funetions, 
and 3) ereative teleseoping algorithms from [ |^ Ch. 3]. Thus, it enables automated deduetion 
of differential equations for holonomie funetions (e.g., our m.g.f. infinite series), their Laplaee 
transform (e.g., our p.d.f.), and their integrals (e.g., our outage probability and ergodie capaeity). 
Conveniently, its symbolie-eomputation abilit>(^ allows for parameters (e.g., Ar, N, Fi, r, ci). 


C. Computer-Algebra-Aided Derivation 

The Mathematica file with HolonomicFunctions eommands that produee the output 
diseussed and employed below ean be downloaded from [ [3^ . Therein, for example, Grobner 
basis eomputation with the eommand Annihilator yields annihilating operators for expres¬ 


sion e from (67). Further, the eommand GreativeTelescoping yields annihilating 


operators for Gn based on its definition as the inner sum in (67), and for Po{z) based on the 


integral in (59). 


^Under “natural boundary” conditions 
'^Inherited from Mathematica. 


30 
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Note that the partieular funetions that enter the differential equations shown below — i.e., 

M^^{s]z), dsM^^{s]z), z); p^^{t]z), dtp^^{t;z), d:,p^^{t]z), dlp^^{t]z)\ d^Po{z), k = 


0 : 4; d’lC{z), A; = 0 : 6 — arise automatieally from (67) by Grdbner basis eomputation and 


ereative teleseoping, and are revealed with the eommand UnderTheStaircase in [36|. 


The steps and outeomes of the proeedure implemented by the eode in |36| are as follows: 


1) Derive SNR m.g.f. differential equations w.r.t. s and based on (67). Then, [361 reveals 
that the funetion veetor 

z) = z) dsM^^{s;z) dM^^{s]z)y 

satisfies the systems of differential equations w.r.t. s and 2 ; 

5^m(s; z) = 0sm(s; z), 92m(s; z) = 0^m(s; z), 


(72) 


with the 3x3 matriees 0^ and 0^ shown only in [36|, due to spaee limitations. 

2) Using results from Step 1, derive p.d.f. differential equations w.r.t. t and x, based on the 


inverse-Laplaee transform. Then, [36| reveals that the funetion veetor 

p{t;z) = {p^,{t;z) dtp^,{t;z) d^p^,{t; z) d^p^,{t; z)y 

satisfies the systems of differential equations w.r.t. t and z 

dtp{t; z) = z), d^p{t; z) = z), 


(73) 


with the 4x4 matriees and shown in |36| . 

3) Using results from Step 2, derive differential equations w.r.t. x for Po{z) and C{z), based 
on their integral relationships from (59) and ( [^ with p.y-^{t; z). Then, [36| reveals that the 
funetion veetors Po(x) = 5x1 with [po(^)]fc = d^Po{z), k = 0 : A, and c{z) = 7x1 with 
[c(x)]fc = d^C{z), k = 0 : Q, satisfy the systems of differential equations 


dzPo{z) = $^Po(2;), d^c{z) = '^zCiz), 


(74) 


where $^ = 5 x 5 and ^2 = 7 x 7 are eompanion matriees [37 p. 109] shown in [36|. 
The above systems of differential equations enable the HGM-based eomputation of the 
SNR p.d.f., outage probability, and ergodie eapaeity, as shown below. 
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VI. Numerical Results 
A. Description of Parameter Settings and Approaches 

For the channel-matrix mean in unit-norm vector a and vector b with the norm in ( [T^ 
are constructed, according to [[^ Eq. (7.29), p. 299], from array response vector^ as 
1 


a = 


b = 


ViVR 

1 


(1 g-J^cos(eR) g-j7r(ArR-l) cos(eR)^T^ 

(1 g-i^cos(0T) 


K 




(75) 


(76) 


’ \l K+l 

assuming uniform linear antenna arrays with interelement spacing of half of the carrier wave¬ 
length. Above, 6 *r and Oj are, respectively, the angles of arrival and departure of the LoS 
component w.r.t. the antenna broadside directions. Unless stated otherwise, we assume ^r = 30° 


and 9-1 equal to the central angle, 9^, of the transmit-side Laplacian power azimuth spectrum [23 
Eq. (4.2)]. Correlation matrix Rj is computed from the AS and 9^ with 


Eqs. (4-3)-(4-5)]. 

Section VI-Bj below shows results for the Stream-1 outage probability for r = 8.2 dB, which 
corresponds to a symbol error probability of 10“^ for QPSK modulation. Thus, the constellation 


size is M = 4, and we show Pq vs. Eb = Eg/ log 2 M = T^/2. On the other hand. Section VI-C 


shows results for the sum rate, i.e., the sum of the ergodic capacities of all streams, in bits 
per channel use (bpcu), vs. AS, K, and 9t. Also shown are simulation results for maximum- 
likelihood detection (ML). 

Unless stated otherwise, presented results have been obtained by running MATLAB R2 012a, 
in its native fixed precision, on a computer with a 3.4-GHz, 64-bit, quad-corep] processor and 8 
GB of memory. Eor the simulation results (in figure legends: Sim.) we have employed, when 
feasible, = 10® samples of n and H for Q, to produce reliable results for Po as low as 
10“®. Then, series results (in legends: Series) have been produced by truncating (67) as in 


[19 Section V.E], i.e., new terms have been added until: 1) their relative change falls below 


10“^®, or 2) n < = 150, as additional terms in (67) lead to numerical divergence because 

the arising large numbers are represented with poor precision. Numerical divergence is indi¬ 
cated in legends with Series*. Outage probability results for full-Rayleigh fading (in legend: 


'See 


0Fig. 


7.3b, p. 296, Eq. (7.20), p. 297] for geometry and derivation details. 


'^Nevertheless, we have run single instances of MATLAB when measuring the computation time (with tic, toe.) 
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Rayleigh, Exp .) have been obtained with expression Pq = > obtained from (|60|) based 


on RemarkFinally, HGM results (in legends: HGM) have been produeed by solving — with the 
MATLAB ode 4 5 function with tolerance levels of 10“^° — the systems of differential equations 


in (74). Then, for the outage probability, the initial condition Po(-So) has been computed accurately 
with and at zq = 0.05692, which arises from d^ for K = —25 dB, A^r = 6, Nj = 4, 


and Rt = Iat^ • Finally, sum rate results have been obtained by adding the ergodic capacities of 
the streams. 

Results are shown for K and AS values relevant to their lognormal distributions for WINNER 


II scenarios AI (indoors office) and C2 (urban macrocell), under LoS propagation [22 Table 5.5]: 
1) averages of these distributions, i.e., for K = 7 dB, and for AS = 51° and 11°, which yield 
low and high antenna element correlation, i.e., |[Rt]i, 2 | = 0.12 and 0.83, respectively; 2) values 
within the range of most likely values [18 Table 1], or 3) random sample^ 


B. Outage Probability Results 

1) Description of Results for K and AS Relevant to Scenario Al, and for Small Nr and Nj: 
Fig. shows results for AS = 51° and K set to values from 0 dB to the upper limit of the 


range expected with 0.99 probability for scenario Al [ 18 Table 1]. Note that the MATLAB series 
truncation diverges for iT = 14 dB and 21 dBp^ whereas HGM and simulation results agree 
at all K. Thus, HGM enables us to investigate the performance degradation likely to occur in 
practice with increasing K for MIMO ZF under full-Rician fading with r = 1. 

Fig. shows results from averaging also over AS and K from their WINNER II lognormal 
distributions for scenario Al. First, simulation has not been attempted due to the long required 
time. (The computation time is explored in more detail below.) Series truncation does not yield 
useful results because of numerical divergence for the larger K values. Only HGM has yielded 
relatively expeditiously a smooth plot whose unshown continuation at sufficiently large Fb has 
revealed the expected diversity ordei[^of A^ = 3 [ 19 Eq. (46)]. 


'^Then, even computing Rt with (2^ Eqs. (4-3)-(4-5)] is time consuming; nevertheless, the employed 2,100 samples of AS 
and K have yielded smooth outage probability plots. 

*"^Our series truncation in Mathematica, with its arbitrary precision, converged also for fT = 14 dB, but required one hour 
vs. a few seconds for HGM; series truncation in Mathematica was not tried for K = 21 dB. 

15n 


The expected diversity order is also noticeable from the plots for TT = 0 dB and 7 dB in Fig 


0 
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Fig. 1. Stream-1 outage probability for A'r = 6, A^t = 4, AS = 51° (i.e., scenario A1 mean), and various values of K, 
including K — 7 dB (i.e., scenario A1 mean). Series results for K — 14, 21 dB do not appear because of numerical divergence. 


Figs. and 1^ depict the same Fb range in order to reveal that: 1) setting AS and K to their 
averages can substantially overestimate performance vs. averaging over AS and K — compare 
the blue dash-dotted plot in Fig. [^with the solid black plot in Fig. 2) making the assumption 
of full-Rayleigh fading instead of full-Rician fading leads to unrealistic performance expectations 
— compare the plots in Fig. 

2) Description of Results for K, AS Relevant to Scenarios Al, C2, and for Increasing Nr, 
Nt: 


Table [^summarizes results of several numerical experiments for K and AS set to their averages 
for scenarios Al and C2, and for the pair (Ar, Nj) set to Na x (6,4), with Na shown in the second 
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Fig. 2. Stream-1 outage probability for A^r = 6, A'r = 4, averaged also over the WINNER II lognormal distributions of K 
and AS for scenario Al. Results corresponding to Rician, Series do not appear because of numerical divergence. 


TABLE I 

Results for K = 7 dB, AS = 51° (i.e., scenario Al) and AS = 11° (C2), and (A^r, Wt) = Na x (6,4). 


AS 

Na 

Tb (dB) 

Po = [ax 10 fox 10 °1 

Series 

Sim. (N, = 10°) 

HGM 

51° (Al) 

1 

[15, 25J 

a = 1.53,fo = 2.15 

1.3 s / 

31 s 

20 s / 

51° (Al) 

2 

[11, 17J 

a = 1.74, fo = 4.26 

1.3 s ^ 

53 s 

20 s / 

51° (Al) 

5 

[6, 

a= 1.39, fo = 6.39 

1.3 s ;( 

520 s 

20 s / 

51° (Al) 

10 

[2, 4.5J 

a = 2.35,fo = 2.45 

1.3 s ;( 

2,300 s 

20 s / 

51° (Al) 

15 

[0, 2J 

a = 1.98,fo= 1.61 

1.3 s ;( 

8,800 s 

20 s / 

51° (Al) 

100 

[-9.2, -8.5] 

a = 2.72,fo = 2.57 

1.3 s X 

estimated : 1.9 x 10° s X 

20 s / 

11° (C2) 

1 

[23, 32] 

a = 1.12, fo = 3.01 

1.3 s / 

31 s 

20 s / 

11° (C2) 

2 

[18.5, 24.5] 

a = 1.43,fo = 3.36 

1.3 s X 

54 s 

20 s / 

11° (C2) 

10 

[5, 7.5] 

a = 2.12,fo = 2.09 

1.3 s X 

2,400 s 

20 s / 


columr|^ The Tb ranges shown in the third eolumn yield Pq in the order of 10 ^-10 ^,as shown 
in the fourth eolumn. The remaining three eolumns show the aetual or estimated eomputation 


'^’Note that Nr does not necessarily have to be much larger than Nj even in massive MIMO 
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time (in seconds), per Fb value. The marks / and / in the ‘Series’ column denote, respectively, 
successful and unsuccessful (i.e., numerical divergence) series computatiorj^ Further, mark 
X in the ‘Sim.’ column indicates infeasible simulation duration. Finally, mark / in the ‘HGM’ 
column indicates successful HGM-based computation. This table demonstrates that, unlike series 
truncation and simulation, HGM enables reliable, accurate, and expeditious ZF assessments for 
realistic K and even large MIMO. 

Fig. [^characterizes ZF performance for iF = 7 dB and AS = 51°, and for the large-MIMO 
setting with Nr = 100 and At = 20. On the one hand, series truncation does not produce useful 
results; on the other hand, HGM results agree with the simulation results, and we have found 
HGM over 30 times fasteiEB 


C. Ergodic Capacity Results 

The ZF ergodic capacity has been computed, for each stream, for Ar = 6, At = 4, = 5°, 

and Fs = 10 dB by: 1) HGM based on ( |74| ) with shown in [361, 2) simulation (also for ML), 


and 3) the infinite series in (67). Results from the series do not appear in the figures because, 
as for the outage probability, its truncation diverges for realistic values of K. 

Fig. [^demonstrates that increasing AS (decreasing antenna correlation) yields increasing ZF 
sum rate and decreasing ML-ZF rate gap. On the other hand. Fig. [^reveals that increasing K 
yields decreasing ZF sum rate and increasing ML-ZF rate gap, for large AS (e.g., 51°). However, 
other (unshown) results indicate that the ML-ZF gap is decreasing for small AS (e.g., 7°) and 
is constant for medium AS (e.g., 12°). 

Finally, Fig. [^ reveals, for AS = 12° and Oc = 5°, a substantial sum rate decrease with 
decreasing |6 *t —Based on Remark [^ because condition Oj = Oc yields worst performance, it 
must also minimize ||hd,i — Hd, 2 i' 2 ,i||- For larger AS, other (unshown) results have revealed more 
moderate rate gain with increasing |6 *t — 0c\. For very large AS (e.g., 51°), the sum rate remains 
unchanged with increasing |6 't — because large AS yields r 2 ,i ~ 0, i.e., ||hd,i — Hd, 2 r 2 ,i|| ~ 
||hd,i||, which is independent of |6 *t — 6*d. Unshown numerical results from the approximating 


*’For (A'r = 6, Ay = 4) numerical convergence is achieved with n = 134, whereas the other {Nk,Nj) pairs yield n = 
rima = 150. Consequently, MATLAB reports about the same computation time (« 1.3 s) for all cases. 

'*When large Ay yields infeasibly-long simulation, HGM results can be validated by checking the diversity order revealed 
by its Po-vs.-Fb plot. E.g., for Ay = 104 and Ay = 100, we have found its slope magnitude to be near the expected N — t>. 
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Fig. 3. Stream-1 outage probability for A'r = 100, A^t = 20, for = 7 dB and AS = 51° (i.e., averages for scenario Al). 
Results corresponding to Series do not appear because of numerical divergence. 


gamma distribution from Remark have revealed it inaeeurate especially for small Nr, Nj, and 
K. On the other hand, we have found that accuracy improves with smaller |6 't — dc\, which 
corroborates Remark |7l 

VII. Summary, Conclusions, and Future Work 

This paper has provided an exact performance analysis and evaluation of MIMO spatial 
multiplexing with ZF, under transmit-correlated full-Rician fading with LoS component of rank 
r = 1. First, we expressed as infinite series the SNR m.g.f. and p.d.f., as well as performance 
measures, e.g., the outage probability and ergodic capacity. However, their numerical convergence 
has been revealed inherently more problematic with increasing K, Nr, and Nj. Therefore, we 
have applied computer algebra to the derived infinite series and deduced satisfied differential 
equations. They have been used for HGM-based computation. Thus, we have expeditiously 
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Arrp=4,Arj^=6;ii' = 7dB;ej^=30°, Be = ex=5°; rs=10dB. 



Fig. 4. ZF sum rate from HGM and simulation vs. AS, for Ar = 6, Ar = 4, A = 7 dB; also, ML sum rate from simulation. 


produced accurate results for the range of realistic values of K and even for large TYr and N^. 
Consequently, we have been able to assess the substantial performance degradation incurred with 
increasing K for ZF when r = 1. Furthermore, HGM has helped reveal that the performance 
averaged over WINNER II AS and K distributions can be much worse than that for average AS 
and K. Finally, we have been able to evaluate the performance for antenna numbers relevant to 
large MIMO reliably and much more expeditiously than by simulation. 

Based on our experience studying MIMO for Rician fading for ZF in this paper and for 
MRC0 in our ongoing work, we expect that performance measure expressions for larger r and 
other transceiver techniques shall entail multiple infinite series in factors proportional to K, 
Ar and At, which shall diverge numerically for realistic values of these parameters. Alternate 

*^Here, MRC refers to the MIMO technique of transmitting and receiving over the dominant channel mode, as discussed in 

g. 
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7Vrp=4,iVj^=6;AS = 52°;6lj^=30°, = 0^=5°; rs=10dB. 



Fig. 5. ZF sum rate from FIGM and simulation vs. K, for A^r = 6, Nr — 4, AS = 52°; also, ML sum rate from simulation. 


computation with the HGM shall require differential equations. Because their by-hand derivation 
from the infinite series shall be intractable, computer algebra shall be indispensable. 


Appendix 


A. Proof of Lemma 


Based on (39) and (42), we can regard [G 2 ].,i = Ar x 1, as a vector of independent complex¬ 
valued Gaussians with variance of 1/2 for the real and imaginary parts, and means 


IE{[G2]i,i} — \/ir2, E{[G2]j,i} — 0, i — 2 : Ar, 


(77) 


which yield 


I[G2]i,iP 

1/2 

l[G2]2,lP |[G2]AfR,lP 

1/2 +■■■+ 1/2 



~ X2(Ar-1)- 


(78) 

(79) 
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Ar^=4,Arj^=6;ii: = 7dB;0j^=3O°, ^0=5°,AS= 12°; rs=10dB. 



Fig. 6. Sum rate vs. 9i for 6c = 5°, when JVr = 6, A'r — A, K = 7 dB, AS = 12°; also, ML sum rate from simulation. 


Now, because T 2 in (44) is upper triangular, we can write the first column of G 2 = U 2 T 2 as 
[ 62 ].,! = [U 2 ],,i[T 2 ]ii. If we set 

[G 2 I ..1 


[T2]i,1 = ||[G2].,i||, [U2].,1 = 


[G2].,iir 


then 


l[U 2 ]l,lP = 


[G 2 ]i ,1 


|[G2]i,iP + |[G2]2,iP H-h |[G2]ArR,l| 


(80) 


(81) 


Finally, using (78), (79), and the independence of [G 2 ]i,i, i = I '■ A^r, one can show that [411 

|[U 2 ]i,i|' ~ B(l,iVR-1,2x2), 


46 


/)i^l-|[U2]i,i|" ~ B(iVR-l,l,2x2). 


The p.d.f. of Pi is then given by [41 j 
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^ „^o "2! 1^7’'i('i“i-i'('r_'V)(«"ir-.dV^• 

'-V- 

=//33 (t);iVR-l,n 2 +l) 


(82) 


where Nr — 1 ,? 7,2 + 1 ) is the p.d.f. of a variable /Ss rsj B(A^r — l,n 2 + 1). Then, the 77 , 1 th 
moment of f3i is 


, p-X2„^2 _^ 

W} = E = E 


(A^r-I) 


ni 


n 2=0 


772 ! 


712=0 


772 ! (772 + iVR) ni 


(83) 


B. Proof of Lemma 

First, let us consider the Nr x A^r matrix G 2 = (G 2 G 2 


ex 


AfR,7VR 


Gd,2, I 


Mr IaTr 


obtained by joining the Nr x Ni matrix G 2 ~ CXiVp, jvi (Gd, 2 ,lAfR ® IatJ froni (39) — whose 


sole nonzero-mean column is [G 2 ].,i — with the NrX N matrix G 2 ~ CXat^^tv (0, Ink ® Iat)- 


Then, paralleling (44), let us consider its QR decomposition, i.e.. 


G 2 — (G 2 G 2 ) — U 2 T 2 , (84) 

with U 2 = Xr X Nr unitary, i.e., U 2 U 2 = U 2 U 2 = Iatr, and T 2 = Xr x Nr upper triangular 


with positive diagonal elements. By partitioning in (84) and using (44), we can write 


T '2 Ti 2 

G2 = (G2 G2) = (U2 U2) 1 ~ I = (U2T2 U2T12 + U2T22), 

0 T 22 


(85) 


where U 2 = Nr x N satisfies U 2 U 2 = 1 x 7 , T 12 = NjX N, and T 22 = N x N is upper triangular 
with positive diagonal elements. 


Hereafter, let us assume that [G 2 ],,i is given, i.e., [U 2 ],,i set as in (80) is given. Then, the 
distribution of 




G 2 U 2 T 2 = (U 2 U 2 )T 2 

= ([U2]..l [U2].,2 ... [U2].,A7, U2)T2 

is invariant to unitary transformations of the columns [U 2 ].,i, Wi = 2 : Ni and the columns of 
U 2 . Thus, we may rewrite 


U2 = (U2 U2)=([U2].,1 UOP), 


(86) 
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where U° = A^r x (TYr — 1) eomprises fixed orthonormal veetors seleeted to form a basis with 


[U 2 ].,i, and P = (iVR — l)x(iVR — 1) is unitary, Haar-distributed [19 See. III.E], not dependent 
on [U 2 ].,i. Using the first row of U° to define 


q'^ = [U°]i,.-P = lx(iVR-l), 


the first row of U 2 from (^) ean be written as 

[U 2 ]i,. = ([U 2 ]i,i q"^) 


Then, based on U 2 U 2 = Iatr and (88), we ean write 
l = ||[U2]i,.f = |[U2]i,iP + ||qf ^ 


|qf = i-|[U2]m 


From ( [87] ) and ( |89| ) we deduee that the veetor 

q 


q 


||q|| V1-|[U2 ]i,i|2 

is uniformly distributed on the unit sphere 
Finally, beeause we ean write 


(87) 


( 88 ) 


(89) 


(90) 


[U2]l,. ^ ([U2]l,l [U2]l,2 ... [U2]l, 


Ni 


test 


- ([U2]i,i qi 


QNi-I QNi 


[U2]l,.) 

. QNR-l)y 


we have that [U 2 ]i, 2 ,..., [U 2 ]i,iv, are the first iVi — 1 elements of q. Thus, we ean write, by also 

[U 2 ] 1 , 2 P + ---+|[U 2 


using ([^, 


A — 


'2jl,7Vil 


1 -I 1 U 2 ], 

ki 

,1 ■ 
2 _ 

2 

f ■ ■ ■ + 

2 

(kll 

p + ■ ■ ■ + Iqn^ 

-1 

P) + (P + • • ■ + 

P) 


Reealling that ||^ is uniformly distributed, we ean deduee that, eonditioned on [G 2 ],,i, i.e., on 
[U 2 ].,i, random variables /94 and /?2 ® 1 — /^4 have the following distributions [41 j: 

/?4 ~ B(iVi-l,iVR-iVi) = B(iVT-2,Ar), 

/)2 = 1 - /?4 ~ B(iVR - A^I, A^I - 1) = B(iV, Nt - 2). 

Beeause the distribution of (32 does not depend on [U 2 ],,i, we also deduee that (32 is independent 

OfftHl-|[U 2 ]l,l |3 
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C. Derivation of Expressions for xi and X 2 

From Remark the normalized vector bn = = A^t x 1 does not depend on K. Defining 



/ 0 0 

\ . 


R = 

\ 0 Rt,A22 

j — iVj X iVx, 

(91) 

r2,i = 

(1 -rl,ir = 

iVj X 1, 

(92) 


we can write pi from (28) and from (41) as: 


/ii = 6* - b^r2,i = b^r2,i = ||b|| bj;;'r2,i, 

h = b^Rb = ||bf b;;'Rb„. 


(93) 

(94) 


Finally, from (10) we have that ||b|p = KN^Nj/ (ii'+l). From (11) we have that [Rp^] ^ ^ 


oc 


{K +1) and RTir 22 'x (Ff +1), i.e., R oc (iF +1), whereas r 2 ,i defined in (23), i.e., r 2 ,i defined 


1 K. These yield: 


xi ^ [R“^]^ J/iip oc 

(95) 

a:2 ^ b^^R-^b oc KN^N^. 

(96) 
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